

*******************************************************
* (1.0) Effets using the true implementation date 
*******************************************************
	
	*For the placebo test, we average effects over the entire post-LLF period 
	
	gen 	years_since_LLFpost 	= 0 if years_since_LLF != . 
	replace years_since_LLFpost 	= 1 if years_since_LLF >= 0 	

	* make interactions parity x years since LLF leaving out parity 0 moms
	forval i = 1/3 { 
			cap drop p`i'_LLFpost
			gen p`i'_LLFpost = 0 if years_since_LLF != . & parity != . 
			replace p`i'_LLFpost = 1 if years_since_LLF >= 0 & parity == `i' 
		}
		
	* make interactions parity x years since LLF x sons leaving out parity 0 moms
	forval i = 1/3 { 			
		cap drop p`i'_LLFpost_son
		gen p`i'_LLFpost_son = p`i'_LLFpost*hasson 	
		}
	
	cap drop post_son
	gen post_son = years_since_LLFpost*hasson 

	order m_year_birth-year  years_since_LLFpre9 years_since_LLFpre8 years_since_LLFpre7 years_since_LLFpre6 ///
		years_since_LLFpre5 years_since_LLFpre4 years_since_LLFpre3 years_since_LLFpre2 years_since_LLFpre1 years_since_LLFpost ///
		pre_9_son pre_8_son pre_7_son pre_6_son pre_5_son pre_4_son pre_3_son pre_2_son post_son ///		 
		p1_LLFpre9 p1_LLFpre8 p1_LLFpre7 p1_LLFpre6 p1_LLFpre5 p1_LLFpre4 p1_LLFpre3 p1_LLFpre2 p1_LLFpre1 p1_LLFpost ///	
		p2_LLFpre9 p2_LLFpre8 p2_LLFpre7 p2_LLFpre6 p2_LLFpre5 p2_LLFpre4 p2_LLFpre3 p2_LLFpre2 p2_LLFpre1 p2_LLFpost ///	
		p3_LLFpre9 p3_LLFpre8 p3_LLFpre7 p3_LLFpre6 p3_LLFpre5 p3_LLFpre4 p3_LLFpre3 p3_LLFpre2 p3_LLFpre1 p3_LLFpost ///
		p1_LLFpre9_son p1_LLFpre8_son p1_LLFpre7_son p1_LLFpre6_son p1_LLFpre5_son p1_LLFpre4_son p1_LLFpre3_son p1_LLFpre2_son p1_LLFpre1_son p1_LLFpost_son ///
		p2_LLFpre9_son p2_LLFpre8_son p2_LLFpre7_son p2_LLFpre6_son p2_LLFpre5_son p2_LLFpre4_son p2_LLFpre3_son p2_LLFpre2_son p2_LLFpre1_son p2_LLFpost_son ///	 
		p3_LLFpre9_son p3_LLFpre8_son p3_LLFpre7_son p3_LLFpre6_son p3_LLFpre5_son p3_LLFpre4_son p3_LLFpre3_son p3_LLFpre2_son p3_LLFpre1_son p3_LLFpost_son 
	
			
	global controls1 	"m_year_birth i.educ age_marriage U5MR_5yr_ave"
	global controls2 	"grain agric_production gdp pop_total_agric"
	
		xi: reg delivered years_since_LLFpre8-years_since_LLFpre2 years_since_LLFpost ///
			pre_8_son-pre_2_son post_son ///
			p1_LLFpre8-p1_LLFpre2 p1_LLFpost  ///
			p2_LLFpre8-p2_LLFpre2 p2_LLFpost ///
			p2_LLFpre8_son-p2_LLFpre2_son p2_LLFpost_son ///
			p3_LLFpre8-p3_LLFpre2 p3_LLFpost ///
			p3_LLFpre8_son-p3_LLFpre2_son p3_LLFpost_son ///
			p_1 p_2 p_3 hasson  p2_son p3_son culturalrev  $controls1 $controls2 i.FE_t i.FE_prov ///
			if year!=1967 & year !=1968 & years_since_LLF<8  & rural==1, cluster(prov)
	keep if e(sample)
	
	*****************************************************
	* (1.2) Grab the Coefficients and CIs
	*****************************************************
	mat coef = e(b)
	mat var = e(V)
		mat coef_all = J(1,7, .) 
		
		local j = 1
		forval i = 2/8 {	
			local k = abs(`i' - 10)
			lincom p`j'_LLFpre`k' + pre_`k'_son 
			local est = r(estimate)
			local se = r(se)
			local dof = r(df)
			local zstat = (`est' - 0) / `se' 
			local pval = 2 * normprob(-abs(`zstat')) 
			local cil = `est' - invttail(`dof', .025)*`se'
			local ciu = `est' + invttail(`dof', .025)*`se'
	
			mat temp1 = `i', `j', 0, coef[1,"p`j'_LLFpre`k'"], var["p`j'_LLFpre`k'","p`j'_LLFpre`k'"], ., .
			mat temp2 = `i', `j', 1, `est', . , `cil', `ciu'
			mat coef_all = coef_all \ temp1 \ temp2
		}
			
		lincom p`j'_LLFpost + post_son 
			local est = r(estimate)
			local se = r(se)
			local dof = r(df)
			local zstat = (`est' - 0) / `se' 
			local pval = 2 * normprob(-abs(`zstat')) 
			local cil = `est' - invttail(`dof', .025)*`se'
			local ciu = `est' + invttail(`dof', .025)*`se'
	
			mat temp1 = 9, `j', 0, coef[1,"p`j'_LLFpost"], var["p`j'_LLFpost","p`j'_LLFpost"], ., .
			mat temp2 = 9, `j', 1, `est', . , `cil', `ciu'
			mat coef_all = coef_all \ temp1 \ temp2
			
		forval j = 2/3 {
		
		forval i = 2/8 {	
			local k = abs(`i' - 10)
			lincom p`j'_LLFpre`k' + pre_`k'_son + p`j'_LLFpre`k'_son
			local est = r(estimate)
			local se = r(se)
			local dof = r(df)
			local zstat = (`est' - 0) / `se' 
			local pval = 2 * normprob(-abs(`zstat')) 
			local cil = `est' - invttail(`dof', .025)*`se'
			local ciu = `est' + invttail(`dof', .025)*`se'
	
			mat temp1 = `i', `j', 0, coef[1,"p`j'_LLFpre`k'"], var["p`j'_LLFpre`k'","p`j'_LLFpre`k'"], ., .
			mat temp2 = `i', `j', 1, `est', . , `cil', `ciu'
			mat coef_all = coef_all \ temp1 \ temp2
			}
			
		lincom p`j'_LLFpost + post_son + p`j'_LLFpost_son
			local est = r(estimate)
			local se = r(se)
			local dof = r(df)
			local zstat = (`est' - 0) / `se' 
			local pval = 2 * normprob(-abs(`zstat')) 
			local cil = `est' - invttail(`dof', .025)*`se'
			local ciu = `est' + invttail(`dof', .025)*`se'
	
			mat temp1 = 9, `j', 0, coef[1,"p`j'_LLFpost"], var["p`j'_LLFpost","p`j'_LLFpost"], ., .
			mat temp2 = 9, `j', 1, `est', . , `cil', `ciu'
			mat coef_all = coef_all \ temp1 \ temp2	
		} 
					
		forval i = 2/8 {	
			local k = abs(`i' - 10)
	
			mat temp = `i', 0, ., coef[1,"years_since_LLFpre`k'"], var["years_since_LLFpre`k'","years_since_LLFpre`k'"], ., .
			mat coef_all = coef_all \ temp
			}	
	
			mat temp = 9, 0, ., coef[1,"years_since_LLFpost"], var["years_since_LLFpost","years_since_LLFpost"], ., .
			mat coef_all = coef_all \ temp
		

*******************************************************
* (2.0) Randomly assign implementation dates 
*******************************************************
	
	forval i = 1/1000 { 
		cap drop temp1
		cap drop temp2

		**************************************************
		* (2.1) Generate placebo treatment years
		**************************************************
		cap drop LLFyear_`i'
		gen temp1 = runiform()
	
		gen LLFyear_`i' = . 
		
		replace LLFyear_`i' = 1964 if temp1 < 1/16 & LLFyear != .
		replace LLFyear_`i' = 1965 if temp1 >= 1/16 & temp1 < 2/16 & LLFyear != .
		replace LLFyear_`i' = 1966 if temp1 >= 2/16 & temp1 < 3/16 & LLFyear != .
		replace LLFyear_`i' = 1967 if temp1 >= 3/16 & temp1 < 4/16 & LLFyear != .
		replace LLFyear_`i' = 1968 if temp1 >= 4/16 & temp1 < 5/16 & LLFyear != .
		replace LLFyear_`i' = 1969 if temp1 >= 5/16 & temp1 < 6/16 & LLFyear != .
		replace LLFyear_`i' = 1970 if temp1 >= 6/16 & temp1 < 7/16 & LLFyear != .
		replace LLFyear_`i' = 1971 if temp1 >= 7/16 & temp1 < 8/16 & LLFyear != .
		replace LLFyear_`i' = 1972 if temp1 >= 8/16 & temp1 < 9/16 & LLFyear != .
		replace LLFyear_`i' = 1973 if temp1 >= 9/16 & temp1 < 10/16 & LLFyear != .
		replace LLFyear_`i' = 1974 if temp1 >= 10/16 & temp1 < 11/16 & LLFyear != .
		replace LLFyear_`i' = 1975 if temp1 >= 11/16 & temp1 < 12/16 & LLFyear != .
		replace LLFyear_`i' = 1976 if temp1 >= 12/16 & temp1 < 13/16 & LLFyear != .
		replace LLFyear_`i' = 1977 if temp1 >= 13/16 & temp1 < 14/16 & LLFyear != .
		replace LLFyear_`i' = 1978 if temp1 >= 14/16 & temp1 < 15/16 & LLFyear != .
		replace LLFyear_`i' = 1979 if temp1 >= 15/16 & LLFyear !=  .
		
		cap drop temp1
		} 
	
		forval iter = 1/1000 { 
		
		**************************************************
		* (2.2) Make the needed interaction terms
		**************************************************
			preserve
			
				keep year LLFyear_`iter' delivered p_1 p_2 p_3 hasson p2_son p3_son culturalrev noson ///
				prov FE_t FE_prov m_year_birth educ age_marriage grain agric_production gdp pop_total_agric U5MR_5yr_ave ///
				parity rural		
	
				cap drop years_since_LLF*
				gen years_since_LLF = year - LLFyear_`iter'
			
				drop if years_since_LLF < -9 | years_since_LLF > 8
				drop if year >=1980
				drop if year <1964
			
				forval i = 1/8 {
					gen 	years_since_LLFpre`i' 	= 0 if years_since_LLF != . 
					replace years_since_LLFpre`i' 	= 1 if years_since_LLF == -`i'
					}

		
				* remake interactions parity x years since LLF leaving out parity 0 moms
				forval i = 1/3 { 
					
					forval j = 1/8 { 
						
						cap drop  p`i'_LLFpre`j'
						gen p`i'_LLFpre`j' = 0 if years_since_LLF != . & parity != . 
						replace p`i'_LLFpre`j' = 1 if years_since_LLF == -`j' & parity == `i' 
						} 
						
						cap drop p`i'_LLFpost
						gen p`i'_LLFpost = 0 if years_since_LLF != . & parity != . 
						replace p`i'_LLFpost = 1 if years_since_LLF >= 0 & parity == `i' 
						
					}
			
				* remake interactions parity x years since LLF x sons leaving out parity 0 moms
				forval i = 1/3 { 
					
					forval j = 1/8 { 
						cap drop p`i'_LLFpre`j'_son
						gen p`i'_LLFpre`j'_son =  p`i'_LLFpre`j'*hasson 
						} 
						
						cap drop p`i'_LLFpost_son
						gen p`i'_LLFpost_son = p`i'_LLFpost*hasson 
						
					}
		
				forval j = 1/8 { 
					cap drop pre_`j'_son
					gen pre_`j'_son = years_since_LLFpre`j' *hasson
					} 
					
					gen 	years_since_LLFpost 	= 0 if years_since_LLF != . 
					replace years_since_LLFpost 	= 1 if years_since_LLF >= 0 	

					cap drop post_son
					gen post_son = years_since_LLFpost*hasson 
						
	
				order m_year_birth-year  years_since_LLFpre8 years_since_LLFpre7 years_since_LLFpre6 ///
					years_since_LLFpre5 years_since_LLFpre4 years_since_LLFpre3 years_since_LLFpre2 years_since_LLFpre1 years_since_LLFpost ///
					pre_8_son pre_7_son pre_6_son pre_5_son pre_4_son pre_3_son pre_2_son post_son ///					 
					p1_LLFpre8 p1_LLFpre7 p1_LLFpre6 p1_LLFpre5 p1_LLFpre4 p1_LLFpre3 p1_LLFpre2 p1_LLFpre1 p1_LLFpost ///					 
					p2_LLFpre8 p2_LLFpre7 p2_LLFpre6 p2_LLFpre5 p2_LLFpre4 p2_LLFpre3 p2_LLFpre2 p2_LLFpre1 p2_LLFpost ///					 
					p3_LLFpre8 p3_LLFpre7 p3_LLFpre6 p3_LLFpre5 p3_LLFpre4 p3_LLFpre3 p3_LLFpre2 p3_LLFpre1 p3_LLFpost ///					 
					p1_LLFpre8_son p1_LLFpre7_son p1_LLFpre6_son p1_LLFpre5_son p1_LLFpre4_son p1_LLFpre3_son p1_LLFpre2_son p1_LLFpre1_son p1_LLFpost_son ///
					p2_LLFpre8_son p2_LLFpre7_son p2_LLFpre6_son p2_LLFpre5_son p2_LLFpre4_son p2_LLFpre3_son p2_LLFpre2_son p2_LLFpre1_son p2_LLFpost_son ///
					p3_LLFpre8_son p3_LLFpre7_son p3_LLFpre6_son p3_LLFpre5_son p3_LLFpre4_son p3_LLFpre3_son p3_LLFpre2_son p3_LLFpre1_son p3_LLFpost_son 
	
		**************************************************
		* (2.3) Estimate placebo regression
		**************************************************

		xi: reg delivered years_since_LLFpre8-years_since_LLFpre2 years_since_LLFpost ///
			pre_8_son-pre_2_son post_son ///
			p1_LLFpre8-p1_LLFpre2 p1_LLFpost  ///
			p2_LLFpre8-p2_LLFpre2 p2_LLFpost ///
			p2_LLFpre8_son-p2_LLFpre2_son p2_LLFpost_son ///
			p3_LLFpre8-p3_LLFpre2 p3_LLFpost ///
			p3_LLFpre8_son-p3_LLFpre2_son p3_LLFpost_son ///
			p_1 p_2 p_3 hasson  p2_son p3_son culturalrev  $controls1 $controls2 i.FE_t i.FE_prov ///
			if year!=1967 & year !=1968 & years_since_LLF<8, cluster(prov)


		**************************************************
		* (2.4) Compile results from placebo regression
		**************************************************
			mat coef = e(b)
				mat coef_`iter' = J(1,1, .) 
				
				local j = 1
				forval i = 2/8 {	
					local k = abs(`i' - 10)
					lincom p`j'_LLFpre`k' + pre_`k'_son 
					local est = r(estimate)
					local se = r(se)
			
					mat temp1 = coef[1,"p`j'_LLFpre`k'"]
					mat temp2 = `est'
					mat coef_`iter' = coef_`iter' \ temp1 \ temp2
				}
					
				lincom p`j'_LLFpost + post_son 
					local est = r(estimate)
					local se = r(se)
			
					mat temp1 = coef[1,"p`j'_LLFpost"]
					mat temp2 = `est'
					mat coef_`iter' = coef_`iter' \ temp1 \ temp2
					
			
				forval j = 2/3 {
				
				forval i = 2/8 {	
					local k = abs(`i' - 10)
					lincom p`j'_LLFpre`k' + pre_`k'_son + p`j'_LLFpre`k'_son
					local est = r(estimate)
					local se = r(se)
			
					mat temp1 = coef[1,"p`j'_LLFpre`k'"]
					mat temp2 = `est'
					mat coef_`iter' = coef_`iter' \ temp1 \ temp2
					}
					
				lincom p`j'_LLFpost + post_son + p`j'_LLFpost_son
					local est = r(estimate)
					local se = r(se)
			
					mat temp1 = coef[1,"p`j'_LLFpost"]
					mat temp2 = `est'
					mat coef_`iter' = coef_`iter' \ temp1 \ temp2
					
				} 
							
				forval i = 2/8 {	
					local k = abs(`i' - 10)
					lincom 0 + years_since_LLFpre`k'
					local est = r(estimate)
					local se = r(se)
			
					mat temp = coef[1,"years_since_LLFpre`k'"]
					mat coef_`iter' = coef_`iter' \ temp
					}
					

				lincom 0 + years_since_LLFpost
					local est = r(estimate)
					local se = r(se)
			
					mat temp = coef[1,"years_since_LLFpost"]
					mat coef_`iter' = coef_`iter' \ temp
					
					
			mat coef_all = coef_all, coef_`iter'
		
		restore
		}

*******************************************************
* (3.0) Plot results
*******************************************************
	preserve
		clear 
		svmat coef_all
		
		rename coef_all1 EventYear
			replace EventYear = EventYear-9
			replace EventYear = EventYear - 1 if EventYear <= -1
		rename coef_all2 Parity
		rename coef_all3 Son
		rename coef_all4 Main
		rename coef_all5 var
		rename coef_all6 cil
		rename coef_all7 ciu
		
			replace cil = Main - invnormal(.025)*sqrt(var) if cil==.
			replace ciu = Main + invnormal(.025)*sqrt(var) if ciu==.
				
		replace Parity = Parity + 1

		forval i = 8/1007 {
			local j = `i'-7
			rename  coef_all`i' coef_`j'
			}	
		
			drop if EventYear < -8 
			drop if EventYear > 7 
			
		order EventYear Parity Son Main var cil ciu coef_*
			
			drop if EventYear < -8 
			drop if EventYear > 8 

		reshape long coef_ , i(EventYear Parity Son Main var cil ciu) j(iteration)
		rename coef_ coef			


		**********************************
		* (3.2) Second order births
		**********************************
		sum Main if EventYear == 0 & Parity == 2 & Son == 1
		local main = r(mean)
		
		cap drop tag
		gen tag = 1 if coef<`main' & EventYear == 0 & Parity == 2 & Son == 1
		count if tag == 1 
		local pval = r(N)/1000
			
		if `pval' == 0 { 
			local pval = "0.000"
			} 
		else if `pval' != 0 { 
			local pval = "0`pval'"
			}
			
		twoway  ///
				(kdensity coef if Parity == 2 & EventYear == 0 & Son == 1) ///			
				, xline(`main', lcolor(red)) xlab(-.08(.02).06) ///
				xtitle("Percentage Point Change in Probability of Birth")  ytitle("") ///
				legend(off) title("Randomly Assigned Treatment Year:" "Second Births with a Male Sibling")	///
				note("P-value from single sided test: `pval'")	

		graph export "$d_fig\Fig_A12a_Permutation.jpg", replace

		sum Main if EventYear == 0 & Parity == 2 & Son == 0
		local main = r(mean)
		
		cap drop tag
		gen tag = 1 if coef<`main' & EventYear == 0 & Parity == 2 & Son == 0
		count if tag == 1 
		local pval = r(N)/1000
		
		if `pval' == 0 { 
			local pval = "0.000"
			} 
		else if `pval' != 0 { 
			local pval = "0`pval'"
			}
			
		twoway  ///
				(kdensity coef if Parity == 2 & EventYear == 0 & Son == 0) ///			
				, xline(`main', lcolor(red))  ///
				xtitle("Percentage Point Change in Probability of Birth")  ytitle("") ///
				legend(off) title("Randomly Assigned Treatment Year:" "Second Births with no Male Sibling")	///
				note("P-value from single sided test: `pval'")	
				
		graph export "$d_fig\Fig_A12b_Permutation.jpg", replace

		**********************************
		* (3.3) Third order births
		**********************************
		sum Main if EventYear == 0 & Parity == 3 & Son == 1
		local main = r(mean)
		
		cap drop tag
		gen tag = 1 if coef<`main' & EventYear == 0 & Parity == 3 & Son == 1
		count if tag == 1 
		local pval = r(N)/1000
			
		if `pval' == 0 { 
			local pval = "0.000"
			} 
		else if `pval' != 0 { 
			local pval = "0`pval'"
			}
			
		twoway  ///
				(kdensity coef if Parity == 3 & EventYear == 0 & Son == 1) ///			
				, xline(`main', lcolor(red)) xlab(-.12(.02).06) ///
				xtitle("Percentage Point Change in Probability of Birth")  ytitle("") ///
				legend(off) title("Randomly Assigned Treatment Year:" "Third Births with a Male Sibling")	///
				note("P-value from single sided test: `pval'")	
				
		graph export "$d_fig\Fig_A12c_Permutation.jpg", replace

		sum Main if EventYear == 0 & Parity == 3 & Son == 0
		local main = r(mean)
		
		cap drop tag
		gen tag = 1 if coef<`main' & EventYear == 0 & Parity == 3 & Son == 0
		count if tag == 1 
		local pval = r(N)/1000
			
		if `pval' == 0 { 
			local pval = "0.000"
			} 
		else if `pval' != 0 { 
			local pval = "0`pval'"
			}
			
		twoway  ///
				(kdensity coef if Parity == 3 & EventYear == 0 & Son == 0) ///			
				, xline(`main', lcolor(red)) xlab(-.08(.02).06) ///
				xtitle("Percentage Point Change in Probability of Birth")  ytitle("") ///
				legend(off) title("Randomly Assigned Treatment Year:" "Third Births with no Male Sibling")	///
				note("P-value from single sided test: `pval'")		
				
		graph export "$d_fig\Fig_A12d_Permutation.jpg", replace

		**********************************
		* (3.4) Fourth order births
		**********************************
		sum Main if EventYear == 0 & Parity == 4 & Son == 1
		local main = r(mean)
		
		cap drop tag
		gen tag = 1 if coef<`main' & EventYear == 0 & Parity == 4 & Son == 1
		count if tag == 1 
		local pval = r(N)/1000
			
		if `pval' == 0 { 
			local pval = "0.000"
			} 
		else if `pval' != 0 { 
			local pval = "0`pval'"
			}
			
		twoway  ///
				(kdensity coef if Parity == 4 & EventYear == 0 & Son == 1) ///			
				, xline(`main', lcolor(red)) xlab(-.14(.02).08) ///
				xtitle("Percentage Point Change in Probability of Birth")  ytitle("") ///
				legend(off) title("Randomly Assigned Treatment Year:" "Fourth Births with a Male Sibling") ///
				note("P-value from single sided test: `pval'")			
				
		graph export "$d_fig\Fig_A12e_Permutation.jpg", replace

		sum Main if EventYear == 0 & Parity == 4 & Son == 0
		local main = r(mean)
		
		cap drop tag
		gen tag = 1 if coef<`main' & EventYear == 0 & Parity == 4 & Son == 0
		count if tag == 1 
		local pval = r(N)/1000
			
		if `pval' == 0 { 
			local pval = "0.000"
			} 
		else if `pval' != 0 { 
			local pval = "0`pval'"
			}
			
		twoway  ///
				(kdensity coef if Parity == 4 & EventYear == 0 & Son == 0) ///			
				, xline(`main', lcolor(red)) xlab(-.12(.02).06) ///
				xtitle("Percentage Point Change in Probability of Birth")  ytitle("") ///
				legend(off) title("Randomly Assigned Treatment Year:" "Fourth Births with no Male Sibling")	///
				note("P-value from single sided test: `pval'")		
				
		graph export "$d_fig\Fig_A12f_Permutation.jpg", replace
